ROC Analysis and a Realistic Model of Heart Rate Variability 
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We have carried out a pilot study on a standard collection of electrocardiograms from patients 
who suffer from congestive heart failure, and subjects without cardiac pathology, using receiver- 
operating-characteristic (ROC) analysis. The scale-dependent wavelet-coefficient standard deviation 
o"wav('7i), a multiresolution-based analysis measure, is found to be superior to two commonly used 
measures of cardiac dysfunction when the two classes of patients cannot be completely separated. A 
jittered integrate-and-fire model with a fractal Gaussian-noise kernel provides a realistic simulation 

QQ ' of heartbeat sequences for both heart-failure patients and normal subjects. 

Q\ \ PACS number(s): 87.10.-he, 87.80.-Hs, 87.90.-hy 

a^ 

Introduction - The interbeat-interval (R-R) time series of the human heart exhibits scahng behavior, as evidenced by 
Ph ' the power-law form of its power spectrum which decreases as f~^ for sufficiently low frequencies / However, other 
^ , features associated with physiological markers |^ are also present in the power spectrum at particular frequencies. 
Moreover, it is well known that the heartbeat time series is nonstationary, reflecting biological adaptability. 

Multiresolution wavelet analysis |4|-H provides an ideal means of decomposing a signal into its components at 
different scales, and at the same time has the salutary effect of ehminating nonstationarities |p|-|lT|. We previously 
carried out a study in which wavelets were used to analyze the sequence of interbeat intervals from a standard 
, electrocardiogram (ECG) database |l3). Using the wavelet-coefficient standard deviation (Twav(w), where m is the 
fSl ' scale, we discovered a critical scale window over which it was possible to perfectly discriminate heart-failure patients 
C^l , from normal subjects. The presence of this scale window has been confirmed in a recent Israeli-Danish study of diabetic 
• patients who had not yet developed clinical signs of cardiovascular disease Q. These two studies, in conjunction 
with earlier investigations involving the counting statistics of the heartbeat (as opposed to the trme-interval statistics 
considered here) [p5|-p7[ , have led us to the conclusion that scale-dependent measures (such as the wavelet-coefficient 
standard deviation) outperform scale-independent ones (such as the scaling exponent S) in discriminating patients 
with cardiac dysfunction from normal subjects. 
^ ' The perfect separation achieved in our initial study endorsed the choice of (Twav(TO) as a measure of importance, at 
least for m = 5 (corresponding to 2^ — 32 heartbeat intervals). The results of most studies are seldom so clear-cut, 
i' however. In circumstances where there is incomplete separation between two classes of subjects, as observed for other 
O , measures using these identical data sets , or in applying our measure to large collections of out-of-sample data 

^ ' sets, the relative abilities of different measures in determining the presence of disease is best established by the use of 
'q ] receiver-operating-characteristic (ROC) analysis [Q. ECG recordings of reduced duration also give rise to incomplete 
• • . separation using our wavelet measure. 
. 5^ ' In this Letter we use ROC analysis to quantitatively compare the tradeoff between data length and discriminability 
, provided by crwav('7i), and by two other widely used heart rate variability measures of cardiac dysfunction. We then 
' develop a mathematical model for heartbeat time-series generation for both heart-failure patients and normal subjects. 
Using ROC analysis to identify cardiac dysfunction - We wish to establish the tradeoff between reduced data length 
on one hand, and misidentifications (misses and false positives) on the other. The ROC curve is a plot of sensitivity vs 
specificity as the threshold parameter is swept; its use requires no assumptions about the statistical nature of the data. 
The area under the ROC curve serves as a well-established index of diagnostic accuracy pO| ; the minimum value (0.5) 
arises from assignment by pure chance whereas the maximum value (1.0) corresponds to perfect assignment. Carrying 
out the ROC calculations permits us to quantitatively compare the abilities of different measures in diagnosing disease. 

In Fig. 1 we present ROC areas, as a function of data length, using the three different measures. The solid 
curve in Fig. 1(a) shows ROC area when discriminability is determined by the wavelet-coefRcient standard deviation 
o'wav(w — 5), using Daubechies 10-tap wavelets Figure 1(b) represents the area when the interbeat-interval 
standard deviation Cint is used instead. The importance of this measure has long been known pit] and it is now 
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commonly used in heart rate variability analysis ||2^. Figure 1(c) provides the ROC area for yet another well-known 
measure, the spectral scaling exponent S estimated at ultra-low (< 0.003 Hz) and at very-low (< 0.04 Hz) frequencies 

ill- 

It is clear from Fig. 1 that crwav(w = 5) is the only measure of the three that, for the ECG recordings in our pilot 
study, ever achieves an ROC area of unity, thereby indicating perfect ability to separate the heart-failure patients 
from the normal subjects, and it does so with as few as 20 000 heartbeats (corresponding to 4 or 5 hours of data). 
It is equally evident from Fig. 1 that the measure of choice for ECG recordings with fewer than 3 500 heartbeats 
(corresponding to about 45 minutes of data) is not crwav('7i — 5) but rather (Tint- This transition occurs because Cint 
depends only on the short-term behavior of the R-R sequence [|l2|,^ whereas the wavelet measure depends on both 
the short- and long-term behavior. It is also apparent from Fig. 1 that (Twav(w = 5) and a-mt always outperform 
the scaling exponent 6, whatever the data length. Moreover, because S reflects the long-duration properties of the 
interbeat-interval sequence its error brackets become unacceptably large as data length decreases (see Fig. 

1(c)) so that it can only be reliably calculated for long data sets. Since we have previously shown ||l^] that the 
spectrum of the R-R sequence, the spectrum of the generalized rate, and the Allan factor all exhibit scaling exponents 
that are similar in value (denoted S, (3, and 7, respectively, in [|l6| ), the use of scale-dependent measures such as 
Cwav(TO = 5) and CTjnt is likely to prove superior to the use of scale-independent measures such as 5, /3, or 7. 

Generating a realistic heartbeat sequence - The generation of a mathematical point process that faithfully emulates 
the human heartbeat could be of importance in a number of venues, including application to pacemaker excitation. 
Integrate-and-fire (IF) models, which are physiologically plausible, have been developed for use in cardiology [§3|,|4| . In 
the paper by Berger et al. p^ , for example, an integrate-and-fire model was constructed by integrating an underlying 
rate function R{t) until it reached a fixed threshold 9, whereupon a point event was triggered and the integrator reset. 
The occurrence time for the (fc 4- l)st beat is then implicitly given hy 9 — J^''^^ dr. Modeling the stochastic 
component of the rate function as band- limited fractal Gaussian noise (FGN), which introduces scaling behavior into 
the heart rate, and setting 9 = 1, results in improved agreement with experiment [|l6|. This fractal-Gaussian-noise 
integrate-and-fire (FGNIF) process requires four parameters: the scaling exponent, the relative strength of the FGN 
spectrum, and lower and upper limits for the noise band. The FGNIF has been quite successful in fitting a whole host 
of interval- and count-based measures of the heartbeat sequence for both heart-failure patients and normal subjects 
p6t . However, it is not able to accommodate the differences observed in the behavior of (Twav(TO) for the two classes 
of data. 

To remedy this defect, we have constructed a jittered version of this model which we dub the FGNJIF. The process 
is generated as follows. Preliminary event occurrence times t^^ are generated by the FGNIF; a Gaussian jitter 
distribution of standard deviation J is then convolved with each of the t^^ to determine the times of the final points 
ti p5[ . Increasing the jitter parameter imparts additional randomness to the R-R time series at small scales, thereby 
increasing dwav at small values of m and, concomitantly, the power spectral density at large values of the frequency 
/• 

In Fig. 2, we present simulations for the wavelet-coefficient standard deviation cr™° versus scale m using the 
FGNJIF model. For J = (circles), the results reduce to those for the FGNIF model, and the shape of the simulated 
curves is in reasonably good accord with experimentally observed curves for normal subjects (Fig. 3(b) left). As 
the jitter standard deviation J increases above 0, the curves in Fig. 2 bend upward at small values of the scale 
m, and the curves begin to match those for heart-failure patients (Fig. 3(b) right). The increased jitter also gives 
rise to a whitening of the spectrum at high frequencies, as expected, so that the distinctions in the spectra between 
heart-failure patients and normal subjects are properly mimicked by the FGNJIF model (compare Fig. 3(d) right 
and left). Typical values of J that accommodate the data lie in the range 0.01 to 0.06 for heart-failure patients and 
in the range to 0.02 for normal subjects. 

The input parameters used to generate some of the simulated curves displayed in Fig. 2, as well as estimates of the 
quantity a obtained from cr^™ over different ranges of m p^ , are presented in the lower portion of Table 1. Care 
should be exercised in referring to a, however, since it is sometimes estimated over a very narrow region of m and 
therefore cannot be considered as a slope or scaling exponent. Examination of Figs. 2 and 3(b) (right) reveals that 
the heart-failure wavelet-coefficient standard deviation curves gently bend upward at small scales. It is clear from 
the lower portion of Table 1 that increasing the jitter at the input to the simulation ( Jinp) results in a reduction 
in a(l < TO < 3), mimicking the observations for heart-failure patients as shown in the upper portion of Table 1. 
A related distinction for the two classes of patients was observed by Peng et al. |l^ using detrended fluctuation 
analysis but, in contrast to us, those authors literally speak of two scaling regions. The effective scaling exponents 
a(3 < m < 10) in the large-scale (low-frequency) regime, slowly decrease with increasing Jinp, as expected, since the 
curves bend up at low scales. 

In Fig. 3 we demonstrate that the FGNJIF simulation does a rather remarkable job of reproducing the actual 
data for a number of key measures used in heart rate variability analysis. The results of the model calculations 
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are compared with ECG data for a single normal subject (left column) and for a single heart-failure patient (right 
column). The input parameters (mean interbeat interval (ri)inp and wavelet scaling exponent a-mp) used in the model 
were obtained from the data sets. The jitter standard deviation used for this particular normal simulation was zero 
whereas it was J = 0.023 for the heart-failure simulation. 

Finally, it is of interest to examine the global performance of the fractal-Gaussian-noise jittered integrate-and-fire 
(FGNJIF) model for the entire collection of data sets in our pilot study. To this end we construct a simulated ROC 
curve using the measure cr^™(r7i — 5). The results for the underlying area are presented as the dashed curve in 
Fig. 1(a). It is derived from 27 simulations |25|, each with 70000 events, to match the pilot-study data. Model 
parameters were determined from the actual ECGs and the scaling exponent was estimated from the slope of the 
wavelet-coefficient standard deviation over the range (3 < m < 10) |2^. The jitter standard deviation value J was 
established by finding the best fit to crwav('Ti)- The corner frequency was taken to be 0.0005 times the mean rate of 
the process (simply 0.0005 for the simulation) as suggested by the results of Rcf. p^ ]. 

The global simulation in Fig. 1(a) (dashed curve) follows the trend of the data (solid curve) quite nicely, but there 
is room for improvement since it falls short of the agreement of the individual simulations illustrated in Fig. 3. It will 
be of interest to consider modifications of the model, including the possibility of nonlinear dynamical behavior p7| , 
that might bring the simulated ROC curves into better accord with the data-based curves. 
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q(1 < m < 10) 


q(1 < m < 3) 


q(3 < m < 10) 


Normal Subjects (12) 








0.79 ±0.08 


1.23 ±0.13 


1.40 ±0.37 


1.22 ± 0.11 


H-F Patients (15) 








0.67 ±0.13 


1.35 ±0.22 


0.26 ±0.60 


1.57 ±0.17 


Interbeat Intervals (Simulations) 


Class 


Input Parameters 


Simulation Results 


('''*) inp 


"inp 


<-^inp 




a(l < m < 10) 


q(1 < m < 3) 


q(3 < m < 10) 


Normal 


0.8 


1.1 


0.00 


0.78 


1.40 


2.03 


1.31 


Heart-Failure 


0.7 


1.4 


0.00 


0.70 


1.67 


2.07 


1.60 








0.01 


0.70 


1.54 


1.28 


1.58 








0.05 


0.70 


1.10 


0.16 


1.35 








0.10 


0.70 


0.84 


0.04 


1.10 



TABLE I. The upper portion of the table presents mean interbeat intervals and quantities a (with their standard deviations) 
estimated from the wavelet-coefficient standard deviation for the 12 normal subjects and 15 heart-failure patients comprising 
our pilot study, collected in two separate groups. Three heart failure patients who also suffer from atrial fibrillation are included 
among the heart-failure patients. Estimates of the effective values of a are determined over the entire range of m, and also 
over two smaller subregions: 1 < m < 3 (the values of a in this region are not scaling exponents - see text) and 3 < m < 10 
[26]. There is a significant difference in the behavior of a for the two classes of patients. The lower portion of the table 
provides estimates obtained from the FGNJIF model. It is clear that the jitter standard deviation J plays a crucial role in 
accommodating the heart-failure observations. 
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FIG. 1. Diagnostic accuracy (area under ROC curve) vs data length (number of heartbeats). A maximum area of unity 
corresponds to the correct assignment of each patient to the appropriate cleiss. The sohd curve in the upper panel (a) is obtained 

by using the wavelet coefficient standard deviation at scale 5 (similar results are obtained at scale 4); the middle panel (b) 
arises from using the interbeat-interval standard deviation; and the lower panel (c) emerges when using the spectral scaling 
exponent. The areas are based on averages of the first 10 data segments for 64, 128, 512, 1024, 3500, and 7000 events (the 
leftmost six data points in (a) and (b)), and on 5, 3, 2, and 1 segments of 14 000, 20000, 35 000, and 70 000 events, respectively 
(the rightmost 4 data points). cTwav is the only measure of the three that achieves 100% sensitivity at 100% specificity for the 
data in our pilot study, and it does so with as few as 20 000 heartbeats (corresponding to 4 or 5 hours of data). For data lengths 
less than 3 500 events (corresponding to about 45 minutes of data), the best performance is provided by Uint rather than (Jwav 
The dashed curve in (a) is derived from 27 simulations of the fractal-Gaussian-noise jittered integrate-and-fire (FGNJIF) model 
(see text). 



FIG. 2. Simulated wavelet-coefficient standard deviation cr^™ curves versus scale m using the FGNJIF model. For J = 0, 
the results reduce to those for the FGNIF model. As the jitter standard deviation J increases from 0, the curves bend up at 
low values of m, accurately mimicking the results for heart-failure patients as is evident from Fig. 3(b) right. 

FIG. 3. Comparison of data from a single normal subject (data set 16265, left column), and a single heart-failure patient 
(data set 6796, right column), with results obtained using the FGNJIF simulation. The parameters (r^) (representing the mean 
interbeat interval) and a (representing the scaling exponent) used in the simulation were obtained from the actual data sets. 
Good fits were obtained by using a jitter parameter J = for the normal subject and J = 0.023 for the heart-failure patient, (a) 
R-R interval sequence over the entire data set. Qualitative agreement is apparent in both cases, (b) Wavelet-coefficient standard 
deviation vs scale. The model reproduces the scaling properties of the data in both cases, particularly the gentle increase of 
CTwav at small scale values, (c) Interbeat-interval histogram. The model captures the narrowing of the histogram (reduction of 
Cint) for heart failure, (d) Spectrum of the sequence of R-R intervals. The simulation captures the subtle distinctions quite 
well, including the whitening of the heart-failure spectrum at high frequencies. 
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Clinical Significance 
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